In [1]:
import pandas as pd
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import osmnx as ox

from numpy import inf
from matplotlib.collections import LineCollection
from mpl_toolkits.basemap import Basemap
%matplotlib inline

In [2]:
def plot_graph(G, m, bbox=None, fig_height=6, fig_width=None, margin=0.02, axis_off=True, equal_aspect=False, bgcolor='w',
               show=True, save=False, close=True, file_format='png', filename='temp', dpi=300, annotate=False,
               node_color='#66ccff', node_size=15, node_alpha=1, node_edgecolor='none', node_zorder=1,
               edge_color='#999999', edge_linewidth=1, edge_alpha=1, use_geom=True):
    """
    Plot a networkx spatial graph.
    Parameters
    ----------
    G : networkx multidigraph
    bbox : tuple
        bounding box as north,south,east,west - if None will calculate from spatial extents of data
    fig_height : int
        matplotlib figure height in inches
    fig_width : int
        matplotlib figure width in inches
    margin : float
        relative margin around the figure
    axis_off : bool
        if True turn off the matplotlib axis
    equal_aspect : bool
        if True set the axis aspect ratio equal
    bgcolor : string
        the background color of the figure and axis
    show : bool
        if True, show the figure
    save : bool
        if True, save the figure as an image file to disk
    close : bool
        close the figure (only if show equals False) to prevent display
    file_format : string
        the format of the file to save (e.g., 'jpg', 'png', 'svg')
    filename : string
        the name of the file if saving
    dpi : int
        the resolution of the image file if saving
    annotate : bool
        if True, annotate the nodes in the figure
    node_color : string
        the color of the nodes
    node_size : int
        the size of the nodes
    node_alpha : float
        the opacity of the nodes
    node_edgecolor : string
        the color of the node's marker's border
    node_zorder : int
        zorder to plot nodes, edges are always 2, so make node_zorder 1 to plot nodes beneath them or 3 to plot nodes atop them
    edge_color : string
        the color of the edges' lines
    edge_linewidth : float
        the width of the edges' lines
    edge_alpha : float
        the opacity of the edges' lines
    use_geom : bool
        if True, use the spatial geometry attribute of the edges to draw geographically accurate edges, rather than just lines straight from node to node
    Returns
    -------
    fig, ax : tuple
    """
    fig, ax = plt.subplots(figsize=(fig_width, fig_height), facecolor=bgcolor)
    node_Xs = [float(node['x']) for node in G.node.values()]
    node_Ys = [float(node['y']) for node in G.node.values()]

    # get north, south, east, west values either from bbox parameter or from the spatial extent of the edges' geometries
    if bbox is None:
        edges = ox.graph_to_gdfs(G, nodes=False, fill_edge_geometry=True)
        west, south, east, north = edges.total_bounds
    else:
        north, south, east, west = bbox

    # if caller did not pass in a fig_width, calculate it proportionately from the fig_height and bounding box aspect ratio
    bbox_aspect_ratio = (north-south)/(east-west)
    if fig_width is None:
        fig_width = fig_height / bbox_aspect_ratio

    # create the figure and axis

    # draw the edges as lines from node to node
    lines = []
    for u, v, key, data in G.edges(keys=True, data=True):
        if 'geometry' in data and use_geom:
            # if it has a geometry attribute (a list of line segments), add them to the list of lines to plot
            xs, ys = data['geometry'].xy
            lines.append(list(zip(xs, ys)))
        else:
            # if it doesn't have a geometry attribute, the edge is a straight line from node to node
            x1 = G.node[u]['x']
            y1 = G.node[u]['y']
            x2 = G.node[v]['x']
            y2 = G.node[v]['y']
            line = [(x1, y1), (x2, y2)]
            lines.append(line)

    # add the lines to the axis as a linecollection
    lc = LineCollection(lines, colors=edge_color, linewidths=edge_linewidth, alpha=edge_alpha, zorder=2)
    ax.add_collection(lc)

    # scatter plot the nodes
    m.scatter(node_Xs, node_Ys, s=node_size, c=node_color, alpha=node_alpha, edgecolor=node_edgecolor, zorder=node_zorder)

    # set the extent of the figure
    margin_ns = (north - south) * margin
    margin_ew = (east - west) * margin


    return fig, ax

In [3]:
def plot_graph_route(G, route, m, bbox=None, fig_height=6, fig_width=None, margin=0.02, bgcolor='w',
                     axis_off=True, show=True, save=False, close=True, file_format='png', filename='temp', dpi=300, annotate=False,
                     node_color='#999999', node_size=15, node_alpha=1, node_edgecolor='none', node_zorder=1,
                     edge_color='#999999', edge_linewidth=1, edge_alpha=1, use_geom=True,
                     origin_point=None, destination_point=None,
                     route_color='r', route_linewidth=4, route_alpha=0.5, orig_dest_node_alpha=0.5,
                     orig_dest_node_size=100, orig_dest_node_color='r', orig_dest_point_color='b'):
    """
    Plot a route along a networkx spatial graph.
    Parameters
    ----------
    G : networkx multidigraph
    route : list
        the route as a list of nodes
    bbox : tuple
        bounding box as north,south,east,west - if None will calculate from spatial extents of data
    fig_height : int
        matplotlib figure height in inches
    fig_width : int
        matplotlib figure width in inches
    margin : float
        relative margin around the figure
    axis_off : bool
        if True turn off the matplotlib axis
    bgcolor : string
        the background color of the figure and axis
    show : bool
        if True, show the figure
    save : bool
        if True, save the figure as an image file to disk
    close : bool
        close the figure (only if show equals False) to prevent display
    file_format : string
        the format of the file to save (e.g., 'jpg', 'png', 'svg')
    filename : string
        the name of the file if saving
    dpi : int
        the resolution of the image file if saving
    annotate : bool
        if True, annotate the nodes in the figure
    node_color : string
        the color of the nodes
    node_size : int
        the size of the nodes
    node_alpha : float
        the opacity of the nodes
    node_edgecolor : string
        the color of the node's marker's border
    node_zorder : int
        zorder to plot nodes, edges are always 2, so make node_zorder 1 to plot nodes beneath them or 3 to plot nodes atop them
    edge_color : string
        the color of the edges' lines
    edge_linewidth : float
        the width of the edges' lines
    edge_alpha : float
        the opacity of the edges' lines
    use_geom : bool
        if True, use the spatial geometry attribute of the edges to draw geographically accurate edges, rather than just lines straight from node to node
    origin_point : tuple
        optional, an origin (lat, lon) point to plot instead of the origin node
    destination_point : tuple
        optional, a destination (lat, lon) point to plot instead of the destination node
    route_color : string
        the color of the route
    route_linewidth : int
        the width of the route line
    route_alpha : float
        the opacity of the route line
    orig_dest_node_alpha : float
        the opacity of the origin and destination nodes
    orig_dest_node_size : int
        the size of the origin and destination nodes
    orig_dest_node_color : string
        the color of the origin and destination nodes
    orig_dest_point_color : string
        the color of the origin and destination points if being plotted instead of nodes
    Returns
    -------
    fig, ax : tuple
    """

    # plot the graph but not the route
    fig, ax = plot_graph(G, m, bbox=bbox, fig_height=fig_height, fig_width=fig_width, margin=margin, axis_off=axis_off, bgcolor=bgcolor,
                         show=False, save=False, close=False, filename=filename, dpi=dpi, annotate=annotate,
                         node_color=node_color, node_size=node_size, node_alpha=node_alpha, node_edgecolor=node_edgecolor, node_zorder=node_zorder,
                         edge_color=edge_color, edge_linewidth=edge_linewidth, edge_alpha=edge_alpha, use_geom=use_geom)

    # the origin and destination nodes are the first and last nodes in the route
    origin_node = route[0]
    destination_node = route[-1]

    if origin_point is None or destination_point is None:
        # if caller didn't pass points, use the first and last node in route as origin/destination
        origin_destination_lats = (G.node[origin_node]['y'], G.node[destination_node]['y'])
        origin_destination_lons = (G.node[origin_node]['x'], G.node[destination_node]['x'])
    else:
        # otherwise, use the passed points as origin/destination
        origin_destination_lats = (origin_point[0], destination_point[0])
        origin_destination_lons = (origin_point[1], destination_point[1])
        orig_dest_node_color = orig_dest_point_color

    # scatter the origin and destination points
    m.scatter(origin_destination_lons, origin_destination_lats, s=orig_dest_node_size,
               c=orig_dest_node_color, alpha=orig_dest_node_alpha, edgecolor=node_edgecolor, zorder=4)

    # plot the route lines
    edge_nodes = list(zip(route[:-1], route[1:]))
    lines = []
    for u, v in edge_nodes:
        # if there are parallel edges, select the shortest in length
        data = min([data for data in G.edge[u][v].values()], key=lambda x: x['length'])

        # if it has a geometry attribute (ie, a list of line segments)
        if 'geometry' in data and use_geom:
            # add them to the list of lines to plot
            xs, ys = data['geometry'].xy
            lines.append(list(zip(xs, ys)))
        else:
            # if it doesn't have a geometry attribute, the edge is a straight line from node to node
            x1 = G.node[u]['x']
            y1 = G.node[u]['y']
            x2 = G.node[v]['x']
            y2 = G.node[v]['y']
            line = [(x1, y1), (x2, y2)]
            lines.append(line)

    # add the lines to the axis as a linecollection
    lc = LineCollection(lines, colors=route_color, linewidths=route_linewidth, alpha=route_alpha, zorder=3)
    ax.add_collection(lc)

    # save and show the figure as specified
    return fig, ax

In [5]:
m = Basemap(resolution='f', urcrnrlat=43.7207,llcrnrlat=43.6255,urcrnrlon=-70.1874,llcrnrlon=-70.3745)
graph_portland = ox.graph_from_place('portland, maine')
fig, ax = plot_graph(graph_portland, m, fig_height=20, fig_width=20, node_alpha=0.0);
m.drawmapboundary(fill_color='cornflowerblue', ax=ax);
m.fillcontinents(color='cornsilk');
fig.savefig('./portland')


---------------------------------------------------------------------------
JSONDecodeError                           Traceback (most recent call last)
/usr/local/lib/python3.5/dist-packages/osmnx/core.py in overpass_request(data, pause_duration, timeout, error_pause_duration)
    282         try:
--> 283             response_json = response.json()
    284             if 'remark' in response_json:

/usr/local/lib/python3.5/dist-packages/requests/models.py in json(self, **kwargs)
    893                     pass
--> 894         return complexjson.loads(self.text, **kwargs)
    895 

/usr/lib/python3.5/json/__init__.py in loads(s, encoding, cls, object_hook, parse_float, parse_int, parse_constant, object_pairs_hook, **kw)
    318             parse_constant is None and object_pairs_hook is None and not kw):
--> 319         return _default_decoder.decode(s)
    320     if cls is None:

/usr/lib/python3.5/json/decoder.py in decode(self, s, _w)
    338         """
--> 339         obj, end = self.raw_decode(s, idx=_w(s, 0).end())
    340         end = _w(s, end).end()

/usr/lib/python3.5/json/decoder.py in raw_decode(self, s, idx)
    356         except StopIteration as err:
--> 357             raise JSONDecodeError("Expecting value", s, err.value) from None
    358         return obj, end

JSONDecodeError: Expecting value: line 1 column 1 (char 0)

During handling of the above exception, another exception occurred:

Exception                                 Traceback (most recent call last)
<ipython-input-5-b00a75efb54b> in <module>()
      1 m = Basemap(resolution='f', urcrnrlat=43.7207,llcrnrlat=43.6255,urcrnrlon=-70.1874,llcrnrlon=-70.3745)
----> 2 graph_portland = ox.graph_from_place('portland, maine')
      3 fig, ax = plot_graph(graph_portland, m, fig_height=20, fig_width=20, node_alpha=0.0);
      4 m.drawmapboundary(fill_color='cornflowerblue', ax=ax);
      5 m.fillcontinents(color='cornsilk');

/usr/local/lib/python3.5/dist-packages/osmnx/core.py in graph_from_place(query, network_type, simplify, retain_all, truncate_by_edge, name, which_result, buffer_dist, timeout, memory, max_query_area_size, clean_periphery)
   1569     # create graph using this polygon(s) geometry
   1570     G = graph_from_polygon(polygon, network_type=network_type, simplify=simplify, retain_all=retain_all,
-> 1571                            truncate_by_edge=truncate_by_edge, name=name, timeout=timeout, memory=memory, max_query_area_size=max_query_area_size, clean_periphery=clean_periphery)
   1572 
   1573     log('graph_from_place() returning graph with {:,} nodes and {:,} edges'.format(len(list(G.nodes())), len(list(G.edges()))))

/usr/local/lib/python3.5/dist-packages/osmnx/core.py in graph_from_polygon(polygon, network_type, simplify, retain_all, truncate_by_edge, name, timeout, memory, max_query_area_size, clean_periphery)
   1473 
   1474         # get the network data from OSM,  create the buffered graph, then truncate it to the buffered polygon
-> 1475         response_jsons = osm_net_download(polygon=polygon_buffered, network_type=network_type, timeout=timeout, memory=memory, max_query_area_size=max_query_area_size)
   1476         G_buffered = create_graph(response_jsons, name=name, retain_all=True, network_type=network_type)
   1477         G_buffered = truncate_graph_polygon(G_buffered, polygon_buffered, retain_all=True, truncate_by_edge=truncate_by_edge)

/usr/local/lib/python3.5/dist-packages/osmnx/core.py in osm_net_download(polygon, north, south, east, west, network_type, timeout, memory, max_query_area_size)
    584             query_template = '[out:json][timeout:{timeout}]{maxsize};(way["highway"]{filters}(poly:"{polygon}");>;);out;'
    585             query_str = query_template.format(polygon=polygon_coord_str, filters=osm_filter, timeout=timeout, maxsize=maxsize)
--> 586             response_json = overpass_request(data={'data':query_str}, timeout=timeout)
    587             response_jsons.append(response_json)
    588         log('Got all network data within polygon from API in {:,} request(s) and {:,.2f} seconds'.format(len(polygon_coord_strs), time.time()-start_time))

/usr/local/lib/python3.5/dist-packages/osmnx/core.py in overpass_request(data, pause_duration, timeout, error_pause_duration)
    298             else:
    299                 log('Server at {} returned status code {} and no JSON data'.format(domain, response.status_code), level=lg.ERROR)
--> 300                 raise Exception('Server returned no JSON data.\n{} {}\n{}'.format(response, response.reason, response.text))
    301 
    302         return response_json

Exception: Server returned no JSON data.
<Response [200]> OK
<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN"
    "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
<head>
  <meta http-equiv="content-type" content="text/html; charset=utf-8" lang="en"/>
  <title>OSM3S Response</title>
</head>
<body>

<p>The data included in this document is from www.openstreetmap.org. The data is made available under ODbL.</p>
<p><strong style="color:#FF0000">Error</strong>: runtime error: open64: 2 No such file or directory /osm3s_v0.7.54_osm_base Dispatcher_Client::1 </p>

</body>
</html>

Driving distance vs Geometric distance


In [5]:
#indexing goodies
nodes = graph_harpswell.nodes()
N = len(nodes)
indices = [i for i in range(0, len(nodes))]
node_index_map = dict(zip(nodes, indices))
index_node_map = dict(zip(indices, nodes))

Pairwise driving distance matrix


In [6]:
#built in method uses an iterative dijkstra's algorithm
#I think this caches the iterative shortest paths search
drive_distances = ox.nx.shortest_path_length(graph_harpswell, weight='length')

#convert dictionary to a numpy matrix
road_distance_matrix = np.zeros(N*N).reshape((N,N))
for i, source in zip(indices, nodes):
    for j, destination in zip(indices, nodes):
        try:
            road_distance_matrix[i][j] = drive_distances[source][destination]
        except:
            road_distance_matrix[i][j] = -inf

True distance matrix


In [7]:
geo_distance_matrix = np.zeros(N*N).reshape((N,N))
for i, source in zip(indices, nodes):
    for j, destination in zip(indices, nodes):
        lat1 = graph_harpswell.node[source]['x']
        lng1 = graph_harpswell.node[source]['y']
        lat2 = graph_harpswell.node[destination]['x']
        lng2 = graph_harpswell.node[destination]['y']
        geo_distance_matrix[i][j] = ox.utils.great_circle_vec(lat1, lng1, lat2, lng2)


/usr/local/lib/python3.5/dist-packages/osmnx/utils.py:267: RuntimeWarning: invalid value encountered in arccos
  arc = np.arccos(cos)

Pandas Dataframe for easy discovery


In [8]:
ids = [(source, destination) for source in nodes for destination in nodes]
road_distance = road_distance_matrix.flatten()
real_distance = geo_distance_matrix.flatten()

df = pd.DataFrame({'id': ids, 'road_dist': road_distance, 'geo_dist': real_distance})
df['difficulty_index'] = df['road_dist']/df['geo_dist']
df = df.fillna(0)
df = df.sort_values('difficulty_index', ascending=False)

mapping


In [58]:
source, destination = df[df['road_dist'] > 1609]['id'].iloc[200]
route = ox.nx.shortest_path(graph_harpswell, source, destination)
m = Basemap(resolution='f', urcrnrlat=43.8840,llcrnrlat=43.7106,urcrnrlon=-69.8847,llcrnrlon=-70.06)
fig, ax = plot_graph_route(graph_harpswell, route, m, fig_height=20, fig_width=20, node_alpha=0.0);
m.drawmapboundary(fill_color='cornflowerblue', ax=ax);
m.fillcontinents(color='cornsilk');
fig.savefig('./harpswell_you_cant_get_there_from_here2')



In [57]:
source, destination = df[df['road_dist'] > 3200]['id'].iloc[300]
route = ox.nx.shortest_path(graph_harpswell, source, destination)
m = Basemap(resolution='f', urcrnrlat=43.8840,llcrnrlat=43.7106,urcrnrlon=-69.8847,llcrnrlon=-70.06)
fig, ax = plot_graph_route(graph_harpswell, route, m, fig_height=20, fig_width=20, node_alpha=0.0);
m.drawmapboundary(fill_color='cornflowerblue', ax=ax);
m.fillcontinents(color='cornsilk');
fig.savefig('./harpswell_you_cant_get_there_from_here')



In [ ]: